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Abstract 

We consider here the moiphogenesis (pattern formation) problem for 
some genetic network models. First, we show that any given spatio-temporal 
pattern can be generated by a genetic network involving a sufficiently large 
number of genes. Moreover, patterning process can be performed by an ef- 
fective algorithm. We also show that Turing's or Meinhardt's type reaction- 
diffusion models can be approximated by genetic networks. 

These results exploit the fundamental fact that the genes form functional 
units and are organised in blocks (modular principle). Due to this modular 
organisation, the genes always are capable to construct any new patterns 
and even any time sequences of new patterns from old patterns. Computer 
simulations illustrate analytical results. 

1 Introduction 

In this paper we consider pattern formation problem in the developmental biol- 
ogy. Mathematical approaches to this problem start with the seminal work by A. 
M. Turing [ 36 ] devoted to pattern formation from a spatially uniform state. Tur- 
ing's model is a system of two reaction-diffusion equations. After [361, similar 
phenomenological models were studied by numerous works (see lfT9ll23l for the 
review). Computer simulations based on this mathematical approach give pat- 
terns similar to really observed ones [ 19 1. However, there is no direct evidence of 
Turing's patterning in any developing organism ([42 1, p. 347). The mathematical 
models are often selected to be mathematically tractable and they do not take into 
account actual experimental genetic information. 
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Moreover, within the framework of the Turing-Meinhardt approach some im- 
portant theoretical questions are left open. For example, whether there exist "uni- 
versal" mathematical models and patterning algorithms that allow to obtain any, 
even very complicated, patterns. In fact, a difficulty in using of simple reaction- 
diffusion models with polynomial or rational nonlinearities is that we have no 
patterning algorithms. To obtain a given pattern, first we choose a reasonable 
model (often using intuitive ideas) and later we adjust coefficients or nonlinear 
terms by numerical experiments (an excellent example of this approach is given 
by the book of H. Meinhardt on pigmentation in shells GUI). 

To overcome this algorithmic difficulty we use genetic circuit models. We 
are going to show that they can serve as "universal models", which are capable 
to generate any spatio-temporal patterns by algorithms. The gene circuits were 
proposed and investigated by many works 1171 fTOl l22l l26l l30l I3T1 l33l l35l Tfor the 
review see Il33l0 in order to use available genetic information, to take into account 
some fundamental properties of gene interaction and understand mechanisms of 
cell gene regulation. 

In this paper we investigate the model from Il2"2"ll2"6l . which is similar to the 
well studied Hopfield neural networks. This model describes activation or depres- 
sion of one gene by another and have the following form: 



where m is the number of genes included in the circuit, jji(x, t) are the concen- 
tration of the i-th protein, Aj are the protein decay rates, Ri are some positive 
coefficients and di are the protein diffusion coefficients. We consider dU) in some 
bounded domain with a boundary d£l. 

The real number measures the influence of the j-th gene on the ?-th one. 
The assumption that gene interactions can be expressed by a single real number 
per pair of genes is a simplification excluding complicated interactions between 
three, four and more genes. Clearly such interactions are possible, however in this 
case the problem becomes mathematically much more complicated. Since the pair 
interaction is capable to produce any patterns, it seems reasonable to restrict our 
consideration only to such interaction. 

The parameters 77^ are activation thresholds and a is a monotone function sat- 
isfying the following assumptions 



dt 



rn 



3=1 



(1) 



a e C°°(R), lim a(z) = 0, lim a(z) = 1, 




(2) 
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I — I < Cexp(-c\z\), a'(0) = l. (3) 
The well known example is a(z) = (1 + tanh z)/2. 

The functions 9i(x) are other activation thresholds depending on x. They can 
be interpreted as densities of proteins associated with the maternal genes. 

This model takes into account only three fundamental processes: (a) decay 
of gene products (the term —XiUi); (b) exchange of gene products between cells 
(the term with A) and (c) gene regulation and protein synthesis. Notice that this 
model of gene circuit can be considered as a Hopfield's neural network [ 14 1 with 
thresholds depending on x and where diffusion is taken into account. The Hopfield 
system is the first model of so-called attractor neural network, both fundamental 
and simple. Analytical methods for the Hopfield models were developed in J5JE21 

CUES). 

Let us fix a function a satisfying ©, © and functions 9%. On the contrary, we 
consider m, K, Aj, di, Ri and rji as parameters to be adjusted. We denote the set of 
these parameters by P : 

P = {m,K,ri,\,d,R}. (4) 

Model ([TJ allows to use data on gene regulation (see ll26ll . where the least 
square approximation of experimental data and simulated annealing were used to 
determine the values of the parameters P). 

In order to study O, many previous works used numerical simulations. For 
example, the work ll3T1 is devoted to the segmentation in Drosophila, in [ 30 ] the 
authors analyse complex patterns occurring under a random choice of the coeffi- 
cients Kij. 

Let us formulate now mathematically our main problem. 

Problem 1.1 (Universal pattern generation problem) Let T > and T < T. 

Given a function z(x,t),x £ 6 [0>^1 an d a positive number t, to find the 
parameters P such that the solution of system with initial conditions yj = 
satisfies 

sup \z(x, t) — yi(x, t)\ < e, xEfl, £ € [T , T]. (5) 

x,t 

Let us consider a biological interpretation of this mathematical formulation. 
We assume that the cell states depend on expression of some genes; we can thus 
identify observed patterns of cell differentiation with gene expression patterns. 
Let us consider these expression patterns as continuous functions z(x,t), where 
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i 6 O, it 6 [0, T] and T > 0. For example, we can assume that z E [0, 1] and if z 
is close to 1, the gene is expressed otherwise the gene is not expressed. 

We consider gene circuits including a single "output" (structural) gene y% = 
Y out and mi "hidden" (regulating) genes. The output gene can change the cell 
states and therefore can predetermine an output pattern z. The hidden genes do 
not influence directly the cell states, they are involved only in an internal cellular 
gene regulation. 

Notice that given pattern z can depend on time t. This fact is important since 
real biological structures are usually dynamical. For stationary patterns z (inde- 
pendent of time) the solution of pattern problem is simple and follows from the 
well known results on neural networks (see Sect. El). 

The main results can be described as follows. 

A) We show that, roughly speaking, any pattern formation process based on 
a reaction-diffusion model can be performed as well by a genetic network, with 
a sufficiently large number of the genes. For each reaction-diffusion model one 
can find an approximating gene network, with the almost same pattern formation 
capacity. This result justifies, to some extent, Turing-Meinhardt's models from a 
genetic point of view. Indeed, these models can be considered as gene circuits. 

B) The second result asserts that, under natural conditions on maternal genes 
densities 0i, the universal pattern generation problem always has a solution. More- 
over, there is a constructive and numerically effective algorithm that allows us to 
find a circuit generating a given pattern. 

Notice that this result is also valid in the absence of diffusion. Indeed, in 
our approach, spatial signalling is not provided by the diffusion process, but by 
space-depending thresholds 6>j. 

Our conditions on maternal gene concentration are necessary and sufficient: if 
they do not hold, it is not possible to approximate any patterns within an arbitrar- 
ily small error. On the contrary, if they are valid, it is possible. If we deal with 
one-dimensional case (for example, we consider a differentiation along anterior- 
posterior axis), then our conditions mean existence of a morphogene gradient 
along this axis. For Drosophila this morphogene is bicoid. So, we show that a 
simple bicoid gradient is capable to produce any chain of complicated time trans- 
formations leading to complex spatial one-dimensional patterns. This result is in 
an agreement with biological observations El S3. To create any two-dimensional 
patterns, we need at least two independent gradients, along anterior-posterior and 
dorso- ventral axes. 

C) We show that the modular organisation and sigmoidal interaction are effec- 
tive tools to form complex hierarchical patterns. 
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Indeed, we show that new, more refined structures, can be obtained by using of 
previous old structures. Also we illustrate that existence of an old structure make 
it easier to produce a new complex one. This property might help to understand 
the usual idea "morphogenesis repeats evolution" ll27ll . see Sect. |4]and |6j 

The paper is organised as follows. In the next section we explain main biologi- 
cal and mathematical ideas beyond these results, in particular, we find connections 
with multilayered network theory and the Hopfield model. In Sect. |3]we describe 
the connection between the reaction-diffusion models and gene circuit systems. 
In Sect. |4]we formulate mathematically pattern formation problem and describe 
main ideas of patterning algorithms. Sect. |5] presents computer simulations il- 
lustrating our analytical results. In particular, as an illustration, we approximate 
numerically, by gene circuit, a reaction-diffusion system for pigmentation of sea 
shells proposed by ll20l . Sect. [6]contains a discussion and concluding remarks. 

All complicated and tedious mathematical details can be found in the Ap- 
pendix. 

2 Main mathematical instruments 

In this section we remind main ideas and results of neural network theory impor- 
tant below. To simplify our statement, we omit some non-essential mathematical 
details (for details, see l38|). 

2.1 Multilayered neural networks 

The neural networks usually consist of a large number of neurons. Each neuron is 
connected to other neurons by directed links with their associated weights. After 
absorbing the inputs, each neuron produces its activation as an output signal to 
other neurons. Each neuron sends a single signal to several neurons at the time. 
Typical problems which may be solved by such nets are pattern classification, 
storing patterns and optimal control problems (see [ 8 ]). 

The simplest example is a single-layer network having one layer of weights. 
The network consists of n input neurons Xj, j = 1, . . . , n and an output neuron 
Y. Each Xj is connected to Y with an associated weight Wj. The output Y is 
given by 




(6) 
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where h is a threshold and a is a strictly monotone function satisfying © and ©. 

Network © can solve only simple classification problems. More powerful, a 
multilayer neural network with p layers consists of one layer of input neurons, an 
output neuron and (p — 1) hidden layers. For p = 2, the corresponding equations 
can be written for instance as 

rn 

Y = a(J2B k z k -h), (7) 

k=l 
n 

z k = <r(y] A kjQj - Vk), (8) 

where qj are states of the input neurons. The remarkable property of this 
network playing the key role in this paper is that any input-output map of the 
form (g 1; q 2 , . . . , q n ) — > 92, • • • , 9n)> where F is a continuous function, can 
be approximated by a network ©-© with a sufficiently large m and appropriate 
weights y4 fcj and S^. 

We shall use below, for brevity, notation 

m 

i=i 

Since a is monotone, the assertion that network ©-(El) approximates any output, 
can be reformulated as follows: by the quantity ^ = Y^k=i B k z k we can approx- 
imate any function, within an arbitrarily small error. This fact results from the 
following well known assertion (see, for example, the works 01E?i n"51l3Sl ). 

Let us consider function ^(q, A, B, 77) of vector argument q = (q x , . . . , q m ) E 
R m depending on the following parameters: the number m > 1, an m x n matrix 
A, and the vectors B and r\ e R m . This function is defined by 

m m 

V{q, A, B, m ,r)) = J2 B k z k = Yl B M A kq ~ Vk)- (10) 

k=l k=l 

We consider this function in a bounded ball consisting of vectors q such 
that \q\ 2 = q{ + . . . + q 2 m < R 2 . 

Lemma 2.1 (Approximation Lemma) If a is monotone and satisfies conditions 
0-0. then for any continuous function Q(q) defined in the ball and for any 



6 



positive number e, there exist a number m > n, matrices A and m-vector r] and 
B such that 

\Q{q)-*S>{q,A,B,m,rj)\<e, qEQ R . (11) 

In other words, given pattern Q and e > 0, we can always find weights A 
and B such that the output of network 0-(0) approximates this pattern, up to 
precision e. 

Approximation (ITU can be obtained by the multilayered network theory (see 
J9j [131 El) or by an application of wavelet type extensions ll38l . For wavelet 
theory see ll2T1l . It is well known that the approximations by dTOb are numerically 
effective as the dimension n of the vector q increases. We know constructive 
algorithms allowing to adjust the parameters A, B, m, 9 (see and references 
therein). 

Notice that the number m of the coefficients B depends polynomially on 
"complexity" of the pattern Q and the precision e. More "complex" the pattern, 
greater m. This complexity can be measured by the integral 01 

Comp(Q) = J \u\\Q(w)\du, (12) 

where Q(cj) is the Fourier transform of the function Q. This means that more 
oscillating functions (patterns) have larger complexities. 

2.2 Large time behaviour of the Hopfield networks 

The Hopfield network [14| is a system of coupled oscillators defined by the dif- 
ferential equations 

d m 

--j- = Ri<r(^2 Kijyj - rji) - y 4 . (13) 
i=i 

Here y,i are neuron states depending on time, Kij is a matrix determining a 
neuron interaction (synaptic matrix), rji are thresholds. Genetic model (UJ) can be 
considered as a generalisation of the Hopfield system such that the neuron states 
and the thresholds depend on a space variable x and the diffusion is taken into 
account. 

We are going to apply some methods developed to investigate attractors of the 
Hopfield neural networks. We recall that dynamics (fl3l is dissipative and thus an 
attractor always exists (on the attractor theory see publications iHl fTTl [Hfl [TH1 12"%1 
34 1 among many others). 
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Dynamics of ( fT3l sharply depends on the synaptic matrix K. If the matrix K 
is symmetric, the attractor usually consists of many equilibria. Such stable large 
time behaviour can be applied to the pattern recognition and associative memory 
problems. 

The large time behaviour of y can become very complex if K is non-symme- 
tric. For instance, depending on K, neuron states can form complicated coherent 
structures that evolve periodically or even chaotically in time. These coherent 
patterns can be described as follows. 

Any (symmetrical or non- symmetrical) m x m matrix K of rank n can be 
represented as a product of two matrices A and B, i.e., 

K = AB, (14) 

where A has size m x n and B has size n x m, 
Let us introduce the new variables 

m 

qi{t) = y %2B lj y j {t) = B l y{t), (15) 

i=i 

where I = 1, 2, . . . , n. 

The dynamical equations for q have the following form 

■^ = -ft + tf { (g,A,B,m,77), (16) 

where ^ are defined by equations similar to dTUb . Time evolution of the new 
variables qi controls the dynamics of all the neuron states y; t . Indeed, we have 

^7T = -yi + cr(Aiq-r}i). (17) 

The functions y(t) can be expressed through q(t) in a simple way by linear 
equations ([T71 ). 

Below we will use new control parameters P, we denote P = {n, m, A, B, 77} 
fixing Ri — 1, \ — 1, di = 0. 

Let us formulate now the following assertion (analogous to the results of ll38l 
39]) describing the complexity of time behaviour of the circuits. 

Lemma 2.2 By the network parameters P, dynamics rti6D can be specified with- 
in an arbitrarily small error. More precisely, for any n, any given continuous 
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functions Qi(q) defined on bounded domain Q, and for any e > 0, we can choose 
parameters P such that 

\Q t + X qi -^i(q,A,B,m,r])\<e, I = 1, 2, . . . , n. (18) 

Therefore, any structurally stable dynamics can be generated by system 

This result shows that the variables qj can exhibit complicated dynamics, pe- 
riodical or chaotical. In particular, any kind of stable chaos can occur in the dy- 
namics of our systems, for example, the Smale horseshoes, Anosov flows, the 
Ruelle-Takens-Newhouse chaos, etc. ll2ll24ll25ll29ll32ll4T1. 

In general, greater the neuron number m, more complex this time dynamics. 
Thus, the neuron states yj also can demonstrate a complicated dynamics however, 
if n « m, all the m neuron states are strongly correlated since they can be 
defined through a relatively small number of the hidden variables. 

For a proof of Lemma 12 .21 see [38 1. 

In the next section we shall show that the gene networks can simulate, in a 
sense, any reaction-diffusion systems. 

Notice that some fundamental and simple biological principles are beyond the 
mathematics. The genes are organised in blocks. The local cell differentiation and 
growth processes are governed by a collective action of these blocks. 



3 Approximation of reaction-diffusion systems by 
gene networks 

We consider, for simplicity, the case of two component reaction-diffusion systems 

Ou 

— = d 1 Au + f(u,v), (19) 
dv 

— = d 2 Av + g(u,v). (20) 

The phenomenological approach based on ([T9l) - (l20l gives excellent results for 
some pattern formation problems (for example such as shell pigmentation |20|), 
where nonlinearities can have the following typical form ll20ll 

u 2 

f = fM{u,v) = av(— - + f3 1 )-Kiu, (21) 

1 + a\U A 
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9 = 9m(u, v) = fc- g + fa) - n 2 v. (22) 

1 + QL\U Z 

We suppose here that all constants a, fy, ai, are positive. 

In these equations, u and v are unknown functions of the space variables x = 
(xi, x 2 , x 3 ) defined in a bounded domain f2. 

System (fT9t- (l20l must be complemented by standard initial and boundary 
conditions. 

Suppose the system of equations that governs patterning is two-component 
system (IT9b — (l20b . where nonlinearities / and g are continuous functions. The 
general multi-component case can be studied in a similar way. Assume solutions 
of (fT9t- (l20l remain globally bounded, i.e., for some positive constants Cj we have 
the estimate 

\u(x,t)\ < Ci, \v(x,t)\ < C 2 , (23) 
for all t > 0, if it holds for t = 0. Let us define the domain Dc lt c 2 as follows: 

Dci,c 2 = {(u,v) : < u < Ci, < v < C 2 }. (24) 

We suppose that initial condition belongs to D Cl: c2 f° r eacn x - 

Our goal is to show that, for a given reaction-diffusion system (fT9t- (l20l we 
can always find an "e- equivalent" circuit ([T]). Namely, for this equivalent circuit 
there exists a smooth map b(y) : (yi, y 2 , ■ ■ ■ , y m ) — > (u, v) transforming the gene 
concentrations to the reagent concentrations and such that time evolution of u, v is 
defined by a new reaction -diffusion system with nonlinearities v), ® 2 (u, v), 
e- close to nonlinearities f(u, v),g(u, v). Roughly speaking we can say that any 
reaction -diffusion system can be realized as a gene circuit. 

To this end, we use Modular Principle. Let us consider a system © having a 
special block structure. Namely, we assume that there exist two kinds of the genes. 
We denote these groups of the genes by y and z, where vector y(x, t) contains m\ 
components and z(x,t) contains m 2 components. Naturally, m = mi + m 2 . We 
consider system (HJ of the special form 

^ = <r(Kf y + Kfz - 9,) + d,A yt , (25) 

3z- — 

= a(K^y + Kfz - 9,) + d 2 Az { . (26) 

Here we use notation © and matrices K. vy , K zz , K zy and K. yz describe inter- 
actions between different groups of the genes. 
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In general, these interactions are not symmetric, i.e., K. yz is not equal to the 
transpose of ~K. zy . 

The coefficients d\ and d 2 coincide with the diffusion coefficients in equations 
<d and <H3). 

We choose the entries of the matrices K. yy , K 22 , K ZJ/ and K. yz as follows: 

#g = *J>i, «S' = 7A, (27) 

and 

K^ = jib jt K% = Uj, (28) 

where a*, Oj, 7^ 7^, fe, are unknown coefficients. 
Let us define "collective variables" 

mi ni2 
U = ^ b iVii V = ^2 b' tZi - ( 29 ^ 

1=1 8=1 

After some calculations ( see the Appendix, part 1) we obtain 

diAu + $i(u,i;), (30) 



c/ 2 Ai; + <l> 2 (u,t;), (31) 

mi 

2j bi(j{aiU + 7if - 6>i), (32) 
i=i 

ni2 

^ b~ia(aiV + - fy). (33) 
i=i 

Applying Lemma 12.11 we notice that for any e > there exist numbers mi, 
777,2, vectors a, 6, a, 6, 7, 7 and 0, such that 

|$i(u,v) — /(u, u)| < e, |$ 2 (u,t;)-0(u,t;)| <e (34) 

for all u, v from some bounded domain. 

This proves the main result of this section: 



du 
~M 

and 

dv 
di 

where 

$2(11, v) ■ 
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Proposition 3.1 Consider problem (Ii9l )-(I20D whose solutions remain in a domain 
Dd,c 2 - Then, if functions f, g are continuous, for any e > 0, there exist such a 
system ([7]) with a sufficiently large number m and coefficients r = (n, r 2 , . . . , r m ) 
and s — (si, s%, . . . , s m ) such that the functions 

m m 

u = ry = r iVii v = -W = Y1 SiVi ^ 

i=l i=l 

satisfy the system 

u t = dxAu + f(u, v) (36) 
v t = d 2 Av + g(u,v), (37) 

where 

\f{u,v)-f(u,v)\<e, (38) 
\g(u,v)-g(u,v)\<e (39) 

for (u,v) G D Cl ,c 2 - 

Therefore, any reaction-diffusion patterning processes on a bounded time in- 
terval [0, T] can be performed as well by genetic networks. In other words, the 
pattern capacity of the gene circuits on bounded time intervals are not less than 
the pattern capacity of reaction-diffusion systems. 

To conclude this section, let us notice that an inverse problem, namely an 
approximation of a neural network by a reaction-diffusion system has been con- 
sidered in @ and OSl . 



4 Programming of spatio-temporal patterns by gene 
circuit models 

In this section we state an analytical algorithm resolving the following problem: 
given spatio-temporal pattern, to find a gene circuit generating this pattern. We 
show that this problem can be solved even without diffusion (di = 0). In our ap- 
proach the space signalling is provided by space-depending activation thresholds. 
It is important from the biological point of view since the molecular transport is 
often performed by non-diffusional mechanisms 1 1 1. For time discrete networks, 
similar results were obtained in ll40ll . 

Beside multilayered network theory (Lemma 12.11) we also use the following 
result. 
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Theorem 4.1 (Superposition Theorem) Let us consider a family T of gene cir- 
cuits with the parameters P 1 , P 2 , . . . , P p , where the functions 9i are fixed and 
identical for all the circuits. Assume these networks generate the output patterns 
Y 1 = y\(x, t),Y 2 = y*(x, t),...,Y' = y*(x, t). 

Then, for any e > and for any continuous positive function F(m, . . . , u p ), 
there is a network generating an output pattern y 1 — Y such that 

\F(Y\x, t), Y 2 (x, *),..., Y p (x, t)) - Y(x, t)\ < e. (40) 

This result can be interpreted as a Superposition Principle. If given circuits 
are capable to produce patterns Y 1 , Y 2 , . . . , Y p , for any function F(u\, . . . , u p ) 
there is a new circuit, which can approximate the pattern z of the form z = 
F(Y 1 , . . . , Y p ), in other words, "superposition by F" of these previous patterns. 
This result also has interesting biological corollaries; we discuss it in Sect. |6] 

Let us describe first the outline of the proof. The proof is based on Modu- 
lar Principle. We suppose that an unknown interaction matrix K of the network 
can be decomposed in blocks. Some blocks contain the known matrices K s cor- 
responding to s-th network of given network family. An additional block deter- 
mines an interaction between new genes and the genes involved in the networks 
of the family T. This structure allows us to apply the approximation results of 
the multilayered network theory Bl|5l[T3l[T5| (see Lemma 12.11) . This assumption 
about the structure of the matrix K also is in agreement with contemporary ideas 
in molecular biology D21EI- The proof (which, by Modular Principle, is quite 
straightforward) can be found in the Appendix. 

Since the basic element of the proof of Superposition Principle is Lemma 12.11 
and the proof of this Lemma gives us an algorithm, therefore we obtain a com- 
plicated but quite constructive algorithm resolving the patterning problem. More- 
over, we can estimate the number of the genes N(z) involved in patterning process 
as a function of the pattern complexity defined by (fT2"h . Namely, using the results 
of the work [3|, we find that N(z) depends polynomially on Comp(z|w), where 
Comply) is a conditional complexity of z respectively given patterns u. To 
explain this relation and its biological meaning, let us consider a simple example. 

Suppose our problem is to construct a periodic one-dimensional pattern z (x) = 
sin kx, where A; is a large number. Our target pattern therefore is sharply oscil- 
lating. Moreover, we have no stored (old) patterns U{ and thus Compel?/) = 
Comp(;z) is proportional to k. In this case, to resolve the pattern approximation 
problem, the network have to involve many genes. 

Assume now that there are old patterns U{ and, in particular, the patterns of the 
form sin k x, cos k x, where k < k but k >> 1. In this case the function z can 
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be expressed through U; L as a polynom of degree P = k/k . Thus Comp(z|w) is 
much less Comp(;z) for large k and k. 

Roughly speaking, a complex target pattern may be simple respectively to 
another complex pattern. We discuss a biological interpretation of this property in 
Sect. |i 

Using Theorem 14.11 we can resolve now the pattern programming problem. 
Suppose the functions 6^(2) possess the following property. They can be consid- 
ered as "coordinates" in the domain fi, i.e., there exist continuous functions gi 
such that 

Xi = 9i(0i(x),6 2 (x), . . . ,0 m (x)), x G Q, i = 1, . . . ,n. (41) 

This condition holds, for example, if m — n and for each i, the function 9i(x) is 
a strictly monotone function of only one variable x, . A biological example can be 
given by the distribution of maternal genes in Drosophila ll42l . 
Let us prove first an auxiliary mathematical result. 

Lemma 4.2 Suppose that condition Mill holds. Then any continuous function 
r \Xi, • • • , X n , t) can be represented as a function ofn + 1 variables 

n = <7(0i(z))(l-exp(-7f)), (42) 

y i = a(^(x))(l-exp(-^)), (43) 

for i — 1 , . . . , n, where k and 7 are two different positive constants. 

Proof: To prove this lemma, let us observe that 

io g y 1 -io g y 1 = /(t), (44) 

where fit) is a strictly monotone function of t. Therefore, t can be written as a 
function of Y\ and Y\. Then any 9i(x) can be presented as a function of Y\ 
and Fi. Using (|4*TT) . one obtains that each X{ IS 3. function of the variables Y s , s = 
1,2, ... ,n and Y x . The lemma is proved. □ 

Let us formulate the main result of this work. This result means that any 
patterning process can be realized by a gene circuit. 

Theorem 4.3 Suppose that condition holds. Then for any continuous posi- 
tive z(x,t),x G Q, t G [0,T], any positive T < T and e there is a system (0) 
such that the solution of this system satisfies the estimate 

\z(x,t) -yi(x,t)\ < e x G Q, t G [T , T]. (45) 



14 



Before start to prove Theorem, let us notice that in the case cfj = (diffusion 
coefficients vanish) condition (HIT) is actually necessary in order to resolve any 
patterning problem. In other words, if it does not hold, there is a pattern, which 
cannot be e- approximated for any e. Indeed, if d{ = 0, solutions of ([TJ are vector 
function y of variables t and 0j. If (|4~T1) does not hold, for some s the pattern 
Ui(x, t) = x s cannot be e-approximated for any e. For d{ ^ condition (HIT) can 
be replaced by a weaker one but we will not consider this question here. 

Proof: 

Theorem I4.3l results from Lemma 14.21 and Theorem 14.11 We take a network 
generating Y\, Y 1 , Y 2 , . . ., Y n . This network has the following structure: 

^ = -72/i+7<r(0i), (46) 

??ji = - K y i + Ka (e i ), (47) 

for i = 1, ...,n. We observe now that ^ = i) and y~i = Yi(x,t). This 
completes the proof. □ 



5 Computer simulations 

We first illustrate the results of Sect. |3l we approximate the reaction-diffusion 
system (TH21)- (l2"2l by a gene network. For this we approximate functions (OTT) and 
d2~2l by sigmoidal functions in order to satisfy inequalities 04l . 

This is a problem of nonlinear approximation since functions $i and $2 de- 
pend linearly on bi and bi but nonlinearly on 7^ 6>j and a^, 7», ^. Coefficients 6j 
and foj can be calculated by the classical least square method, but other coefficients 
have to be determined in a proper way. For instance, we could choose these coef- 
ficients randomly and select the best values (or a satisfying value), however this is 
usually too long when the search space is large. Here this random method cannot 
be used due to the fact that the functions <f2Tb and d22l depend on two variables. 
If we approximate a function of a single variable, this simple random method can 
be useful and we apply it below. 

To make this nonlinear approximation, we use an iterative approach proposed 
by Jones [IV] (see also Barron 0). Jones' result is an iterative version of Approx- 
imation Lemma 12.11 under the conditions of this lemma one can find a sequence 
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of approximating functions (%(q, A, B, r], m)) m€N (denoted ( v I/ m ) mgN for brevity) 
satisfying 

\Q(q)-^ m \ 2 <C/m, (48) 

where C is a constant depending on Q. As already mentioned in Sect. |3 Barron 
[3 1 has related C to the Fourier transform of Q. For more oscillating functions 
Q, the constant C is greater. This constant can be considered as a measure of 
complexity of Q. 

Jones' sequence is defined as follows: \&i = B\cr(A\q — 771), where B\, 
A\, 771 give an almost minimal value of \Q(q) — Then, for any m > 2, 
^ m = a m ^/ m -i + B m cr(A m q - rj m ), where a m , B m , A m and r] m give an al- 
most minimal value of \Q(q) — *& m \. Barron formulates precise conditions 
on these almost minimal values in order to obtain equation (041. but we do not 
use these conditions here. The important point we use is that optimising only one 
sigmoidal function each time, Jones' sequence is able to achieve 0(1 /m) approxi- 
mation. This permits to avoid a global optimisation of all the coefficients involved 
nonlinearly. The linearly involved coefficients a m and B m can be computed by 
the least square method and the nonlinearly involved coefficients A m (which are 
two-dimensional like q) and r] m can be determined by a random method. Such 
approach allows us to approximate functions (I2T1) and (l22l . The numerical pa- 
rameters were a — 8, a% — 1, f3\ — 0, K\ = 2, (3 2 = 1> ^2 = 0. In this case system 
(IT^l-d^t with Neumann boundary conditions has an homogeneous equilibrium 
solution u = (3 2 /ki,v = (^+Q!i/3|)/(q;/3 2 ). On the segment [0,L] withL = 60, 
and with d u — 1, d v — 50, this equilibrium is unstable with respect to some non- 
homogeneous perturbations. Using an initial perturbation on u at the left side 
(u(t — 0, x) = 2uo for x E [0, L/10]), one obtains a non-homogeneous station- 
ary solution (so-called Turing structure). Moreover, the perturbation spreads to 
the right like a wave. See [ 20 1 p. 30. This behaviour is presented in fig. [T] We 
approximated Meinhardt's model by a gene network, with 600 genes (300 to ap- 
proximate function (|2TT) and 300 for (1221 ). The behaviour is qualitatively similar 
to the solution of Meinhardt's model. See fig. |2j 

The second point we illustrate is the universal pattern generation problem: 
given a spatio-temporal pattern, to find a gene network generating pattern. 

We first generate spatio-temporal patterns with one spatial dimension. The 
corresponding gene network can be defined as follows: 

^l = K (a(9 1 (x))-y 1 ), (49) 
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^l = 2n(a(9 1 (x))-y 1 ), (50) 

du ■ — 

= X(R j a(K j y 1 + K$j x - rjj) - uj), j = 1, . . . ,m, (51) 

d m 

-^f- = H°(%2 u j) -2/out). (52) 

3=1 

Notice that diffusion is absent and positional information is provided by the spa- 
ce-dependent threshold 6\ (x) . In this one-dimensional case, the condition on 9i 
means that this function is strictly monotone. 

As it is shown in Sect. |4j t and x are functions of (yi, yi). Namely, 

t = --log(^-l) (53) 



output of Meinhardt's model 




Figure 1: Solution of Meinhardt's model. 
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and 

x = 6- l \a-\—^-)). (54) 

Thus any z{t, x) can be presented as a function of (yi, yi): 

x) = z(time(yi, yi), space(yi, yi)) = Z(yi, yi). (55) 

Notice that y out approximates r = cr(X]j=i Wj) as A — > oo. In turn, r approx- 
imates cr(Y^ =1 Rj &{Tj]ji + Tji/i — 6j)). Hence, to solve the pattern gener- 
ation problem, we have to determine the coefficients Rj, Tj, Tj, 9j such that 
J^jLiRj^iTjUi + T j y 1 — 9j) approximates a~ 1 (Z(y 1 ,y 1 )). It is possible if Z 
is a continuous function. 

The problem is thus to approximate a function of two variables (yi and yi). 
This problem is intractable with the least square and the random methods, and we 
use here again Jones' iterative approximation method (see above). 

To avoid singularities at the lines y x = y x and 2y x = y l5 we have approximated 
this function Z in the image of the bounded rectangle [T , Ti] x [x , xi] by the map 



output of gene network 



1.2 r 




Figure 2: Approximation of Meinhardt's model by a gene network. 
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(t,x) i-> (cr(6 l i(a;))(l - e Kt ), er(0i(:c))(l - e 2k *)), which is a one-to-one map of 
(t,x). 

Fig. |3]and|4]present the output of system (l49l)- (l5TT) approximating the function 

0.1(sin(8t)+sin(16t)) and 0.025(l+tanh(10t-0.5)) sin(8x), respectively, fort £ 
[0, 1] and x £ [0, 1]. We have used 1000 sigmoidal functions for these simulations. 



output of the gene network 
pattern to be generated 




Figure 3: Generation of 0.1(sin(8t) + sin(16£)) by a gene network. 

Also we have generated spatio-temporal patterns with 2 space dimensions. 
The corresponding gene circuit is 

^ = K (a{9 1 {x))-y 1 ), (56) 

dm 

dt 
dyi 
dt 



-K(a(6 2 (x))-y 2 ), (57) 
2 K (a(9 1 (x))-y 1 ), (58) 
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^ = X(R,a(Kj yi + K]y 2 + K 3 y x - Vj ) - Uj ) ,j = l,...,m, (59) 

a m 

^ = X(<r(£u j )-y aat ). (60) 
j'=i 

The time i, the spatial coordinates xi and x 2 can be expressed as functions of 

(2/1,2/2,^1): 

t = --log(^-l), (61) 

K 2/1 

-i = ^(^ 1 (^-),^ 1 (/^-)) (62) 
22/1 - 2/i 2j/i - yi 

and 

/ -1/ 2/i x -1, 2/12/2 xx .™ 

X 2 =g2{(T (tt -),Cr — . (63) 

22/i - 2/i 2yi - 



output of the gene network 
pattern to be generated 




Figure 4: Generation of 0.025(1 + tanh(10t — 0.5)) sin(8:r) by a gene network. 
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Hence, any continuous function z(t, Xi,x 2 ) can be represented as a function of 
{yi, V21 yi)> which has to be approximated by Jones' method in order to solve the 
pattern generation problem. Since t, x\ and x 2 are singular in y\ = %j\ and 2y\ = 
yi, these functions were approximated in the image of the cubic domain [T , T\] x 
[xi,o, a?i,i] x [x 2 ,o, x 2 ,i] by the map (t, x u x 2 ) h-> (cr(0i(:r))(l-e _re *), a(6 2 (x))(l- 
e- Kt ),a(6 1 (x))(l -e- 2Kt )). 

Fig. |5] presents the output of system (I56T)-(|6TH approximating the function 
0.01((xi - 0.5) 2 - (x 2 - 0.5) 2 ) for (xi,ar 2 ) G [0, l] 2 . This function is indepen- 
dent of time, but time-dependent functions have also been approximated (it is not 
shown). We have used 1000 sigmoidal functions for this simulation. 



output of the gene network 
pattern to be generated 




Figure 5: Generation of 0.01((ari ~ 0.5) 2 — (x 2 — 0.5) 2 ) by a gene network. 

The last point we illustrate is the superposition principle and its relation with 
the conditional complexity (see Sect. |4j). The superposition Theorem 14. II states 
that a given network generating a pattern u(t, x) and a given continuous func- 
tion F, one can device a new network generating F(u)(t,x). The number of 
the genes involved in this new network depends on the complexity of the tar- 
get pattern. This complexity can be defined by the Fourier transform of the pat- 
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tern Q. We define the conditional complexity Comp(F(«)(i, x)\u(t, x)) as the 
complexity of F(u)(t,x) considered as a function ofu(t,x). The point is that 
Comp(F(u)(t,x)\u(t,x)) can be much less than Comp(F(u)(t,x)). So gener- 
ating F(u)(t,x) through u(t,x) we may use much less genes than generating 
F(u)(t,x) directly (or, if the same gene number is involved, a better precision 
may be achieved). 

We illustrate this fact by generating cos(8i) for t G [0, 2ir). We produce this 
time function directly and, moreover, we first generate cos(t), then 2 cos 2 (t) — 1 = 
cos(2t), later 2 cos 2 (2t) - 1 = cos(4t) and finally 2 cos 2 (At) - 1 = cos(8t). The 
network generating cos(8t) directly is 



dyi 
dt 



(64) 



du ' 

— = \(Rja(Kjyi - rji) - Uj), j = 1, . . . , m, (65) 



dt 



H^2 u j-Vout), (66) 



where Rj, Kj and r/j are chosen so that J2T=i Rjv(Kjt—rij) approximates cos(8£). 
The network producing cos(8t) indirectly is 

dm _ 1 

dt ' 



9uJ 
dt 



(67) 

X(R]<r(Kj yi //;) '/])../ I mi, (68) 

f = A£>}-y,), (69) 



J'=l 



<9w 2 



dt 



' - \(R)a(K]y % - V f) - u% j = 1, . . . , m 2 , (70) 



= \{R]a{K 2 jyz - ttJ) - j = 1, . . . , m 2 , (72) 
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r) 1112 

^ = A(E«J-*). (73) 



3=1 

= \{R)a(K 2 3 y, - r§) - u% j = 1, . . . , m 2 , (74) 



, j/ =A(£t»}-y OHt ) I (75) 

3=1 

where i^j and r/j have been chosen so that 5^/=a Rj a (Kjt ~ Vj) approxi- 
mates cos(t) and i£| and r/| have been chosen so that £]j=i R 2 jV{K?jX — 77J) 
approximates 2a; 2 — 1 . 

We have used the same number of equations in the two cases, namely 52, (so, 
m = 50 in equation (1651) ) and we have compared the precision achieved. For the 
indirect approximation, we have chose mi = 32 and hence, mi = 5. 

The target pattern is defined by a function of a single variable. In this case 
and with a small number of the genes, using of Jones' method is not obligatory. 
Actually, here by the least square method for the linear coefficients and a random 
choice for the nonlinear ones we achieve a better precision. 

Fig. [6]presents the results. The patterns computed by systems (l64l) - (l66l) and 
(l67ll — (f75l) are denoted respectively "direct approximation" and "indirect approxi- 
mation". 



6 Conclusion 

It is shown that the genetic networks with binary interaction of the genes have a 
formidable patterning capacity. They can produce any spatio-temporal patterns. 
Moreover, it is proved that any reaction-diffusion systems can be approximated 
by genetic circuits. This result allows to connect earlier phenomenological math- 
ematical reaction-diffusion models and more biologically realistic genetic circuits. 

Let us emphasise that, by these circuits, pattern programming can be per- 
formed. This means that, for a given pattern, a circuit that builds this pattern, 
can be found by effective and universal algorithms. One of the most astonishing 
biological revelations of the past twenty years is that much of the basic machinery 
of development is essentially the same, not in all vertebrates but in all the major 
phyla of invertebrates too lH2l . We show therefore that this machinery can be 
described by simple gene circuit models. 
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This pattern programming holds on a basic biological principle: on modular 
organisation of genes. Genes are organised in blocs. Notice that the modular 
principle is confirmed by experimental data of molecular biology (see HI [12J |43 
and references therein). 

We have demonstrated that this modular structure entails an interesting prop- 
erty, which can be named "superposition principle". This superposition property 
means that new patterns can always be obtained by previous (old) patterns. 

As an elementary example explaining a biological interpretation of superpo- 
sition principle we can consider flappers, wings and legs of tetrapodes. It is well 
known that they consist of the same basic elements (numerus, cubitus, radius, 
carpe) but jointed in different ways (see ||2"71 ). Different joinings give wings for 
birds, legs for dogs, flappers for whales etc. When mammals penetrated in water, 
evolution did not invented flappers from zero. Evolution used earlier created pat- 
terns to obtain flappers. Thus, the superposition principle allows us to understand 
why the gene number grows relatively slow in evolution (remind that Drosophila 
has the 14000 genes, C. elegans has the 19000 and Homo sapiens has the 30000 



time 



cos(8t) 
direct approximation 
indirect approximation 



Figure 6: Improvement of approximation by superposition principle. 
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genes). Indeed, as we have explained in Sect. H|and |5J the gene number growth is 
not directly proportional to the pattern complexity; this number is proportional to 
conditional pattern complexity relatively already stored patterns. This conditional 
quantity may be small even if the target pattern is very complex. 

We can thus conclude that the modular organisation of gene interaction leads 
to a minimisation of time and genes in a process of invention of new biological 
structures. A famous basic evolutionary law asserts that the ontogenesis sum- 
marises the philogenesis [27 1. The mathematical results of this paper suggest that 
this law is a direct consequence of gene network organisation. 

A Appendix 

A.l Derivation of equations (l30l) and (I3TI) 

Equations (l25l and (l26t can be rewritten as 

—± = (T (a i u + jiV-9 i ) + d 1 Ay i , (76) 

dz' — 

— f = a(diV + jiU-6i) + d 2 Azi. (77) 
at 

Multiplying the i-th equation in (TToT) by hi and taking the sum over i, we obtain 
CSB and (El). 

A.2 Boundedness of solutions of the Meinhardt equations 

Solutions of (fT9t-(l20b stay bounded: they lie in the domain D Clt c 2 f° r ai l times if 
the corresponding initial data are in this domain. 

Let us show that condition d23t holds with appropriate Ci, C 2 . Let us choose 
such Cj that 

aC 2 {Cl{\ + a^lY 1 + ft) < k x C x (78) 

and 

(3 2 < (aft + n 2 )C 2 . (79) 

First we choose a large C 2 to satisfy (1791) and then we can take a constant 
Ci large enough to satisfy dTHl . Now we can prove that the domain D Cl ,c 2 1S an 
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invariant rectangle for d"HJt-(l2"0"l). This means that on the boundaries u 
v = C 2 the vector field (/, g) is directed inside D Cl ,c 2 - 
This assertion follows from (l78l) and (l79l) . 



= G\ and 



A.3 Proof of Superposition Theorem 

We consider a large network involving a number of genes. First, it involves the 
genes yf, i = 1, . . . , m s , where s = 1, . . . ,p, participating in given networks. 
The corresponding dynamics is defined by the equations 

^l^ = R^(f^K^-9 l (x)-r ] t)-X s l y s l , xefl, t>0, (80) 
i=i 

where s = 1, . . . ,p. 

Moreover, the large network includes additional genes Vk- The time evolution 
of the corresponding concentrations Vk(x, t) is defined by the following equations 

A.' = b k a(J2 M kj y{(x, t) - r) k ) - Xv k , xefl,t>0. (81) 

m 3=1 

At last, the genes v k (x, t) determine the time evolution of the output gene yi 
as follows: 

= Ra (^ S k v k {x, t)) - Ay b xeQ, t>0. (82) 

3=1 

We set the zero initial conditions for all the concentrations 

v k (x, 0) = y*(x, 0) = Vl (x, 0) = 0. 

Let us prove the auxiliary lemma. 

Lemma A.l Given a function z{t) G C^fOjT], positive numbers e and 5 < T, 
there are a function w G C[0, T] and a positive coefficient A such that the solution 
of the Cauchy problem 

( ^jQ = -\X(t)+w(t), X(0)=0, te[0,T] (83) 
satisfies the following inequality 

\X(t)-z(t)\<e, te[8,T\. (84) 
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Proof: The proof of this lemma is elementary. Indeed, let us set 



dz(t) 



+ Xz(t), X = z + X. 



w = 



dt 



Then ( 1531 entails 



dX(t) 



XX (t), X{0) 



z(0). 



dt 



Thus, 

\X(t)\ < \z(0)\exp(-Xt). 



To complete the proof of the lemma, we set 



A>-r 1 log(|^(0)|- 1 e). 



□ 



Notice that the lemma also holds for z £ C[0, T], since any continuous func- 
tion can be approximated by a smooth function. Moreover, if given z is a su- 
perposition of the form z = z(y 1 (t), . . . , y p (t)), where y s are defined by some 
system of autonomous differential equations, then w can also be represented as a 
superposition: w = w{y 1 {t), . . . , y p (t)). 

To finish the proof of Theorem 14.11 it is sufficient to prove that for any contin- 
uous function of the form w(x, t) = w(y\(x, t), . . . , y p (x, £)), where x E f2, t 6 
[0, T] and e > 0, there exists such a choice of the parameters My, A, bk, Sk 
and R in (I8TT) and (1821) that the solutions Vk(x, t) of <f8TT> satisfy the estimate 



for any x £ f2, t £ [^^1- Using the monotonicity of a and choosing a 
sufficiently large R, we simplify the last estimate and obtain 



\W(y{(x,t),...,y P (x,t))-Y,SkV k \<e, x £ Q, t £ [S,T], (86) 




(85) 



fc=i 



fe=i 



where W(x, t) is given. 
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Let us take sufficiently large A > 0. Using ( ISBT ) and (ISTT ). we obtain 

v 

\v k (x,t) - \- x b k a(^M kj y{{x,t) - %)| < e/4. (87) 
j'=i 

Denote = X^S^bk- Now, to finish the proof, it is sufficient to find the 
parameters 

M kj ,7] k ,p k such that 

mo P 

\W(yl(x,t), . . . ,&{x,t)) - X> CT (E M wi(x,t) ~Vk)\< e/4, (88) 

fc=l 3=1 

where x G fi, £ G [5, T\. The existence of this approximation follows from the 

multilayered network theory (see Lemma 12 .lb . 
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